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ABSTRACT 

We solve the spherically symmetric time dependent relativistic Euler equations on a 
Schwarzschild background space-time for a perfect fluid, where the perfect fluid models 
the dark matter and the space-time background is that of a non-rotating supermassive 
black hole. We consider the fluid obeys an ideal gas equation of state as a simple model 
of dark matter with pressure. Assuming out of equilibrium initial conditions we search 
for late-time attractor type of solutions, which we found to show a constant accretion 
rate for the non-zero pressure case, that is, the pressure itself suffices to produce 
stationary accretion regimes. We then analyze the resulting density profile of such 
late-time solutions with the function A/r K . For different values of the adiabatic index 
we find different slopes of the density profile, and we study such profile in two regions: 
a region one near the black hole, located from the horizon up to 50Af and a region 
two from ~ 800M up to ~ 1500M, which for a black hole of 10 9 M Q corresponds to 
<~ O.lpc. The profile depends on the adiabatic index or equivalently on the pressure of 
the fluid and our findings are as follows: in the near region the density profile shows 
values of k < 1.5 and in the limit of the pressure-less case K —> 1.5; on the other hand, 
in region two, the value of n < 0.3 in all the cases we studied. If these results are to be 
applied to the dark matter problem, the conclusion is that, in the limit of pressure-less 
gas the density profile is cuspy only near the black hole and approaches a non-cuspy 
profile at bigger scales within lpc. These results show on the one hand that pressure 
suffices to provide flat density profiles of dark matter and on the other hand show that 
the presence of a central black hole does not distort the density profile of dark matter 
at scales of O.lpc. 
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1 INTRODUCTION 

One of the most important issues related to the dark mat- 
ter problem is that of the dark matter density profiles in 
galaxies. Most of the analyses involve the study of profiles 
based on simulations of structure formation and collapse of 
cold dark matter (CDM). Among the most studied models 
there are for instance the Navarro Frenk White (NFW) den- 
sity profile (|Navarro et a l. 1996, Nav arro et al. 1997|) and 
Moore's model (|Moore et al. 1999|) . Within the study of 
density profiles it is particularly interesting to understand 
the distribution of matter in the innermost region of the 
halo, where different models provide different density slopes 
of the type ~ l/r K ; in particular the NFW and Moore pro- 
files show different behaviors k = 1 and re = 1.5 respectively. 
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Various studies indicate different slopes, for instance in 
( |Klypin et al. 2001| ) it was found a slope with re = 1.5 based 
on the study of halo simulations, from phase space density 
arguments, in ( Taylor fc Navarro 2001[ ) the density profile 
was found to show re = 0.75 instead of re = 1 and resembles 
the NFW profile in the outer region, in (|Colm et al. 2004]) 
based on studies of low mass halos the result was k = 1 
which corresponds to the same limit as NFW, the conclusion 
in (|Diemand et al. 2005) is that CDM halos have cusps with 
re = 1.2, whereas in (|Navarro et al. 2004|) fitting functions 
with re = 0.7 were used for radii or the order of r ~ 0.01 kpc; 
the work in ( Stohcr 2006 ) indicates that within r ~ lkpc the 
solpe corresponds to re ~ 1 and a further extrapolation im- 
plies re ~ which would correspond to a pseudo-isothermal 
model. Also in (|Navarro c t al. 2008J non single slope limits 
are proposed depending on the radial scale considered, for 
instance for r ~ O.lkpc re ~ 0.85, and for r ~ lkpc re = 1.4.. 

Observations on the other hand suggest different pro- 
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files depending on the type of galaxies considered, for in- 
stance dwarf and LSB galaxies are better described by a con- 
stant density core model (jBurkert 19951 IWalter et al. 2008|) . 
specifically, the mass density profile of dwarf galaxies shows 
averages of the order k ~ 0.29 (|Oh et al. 2010[) . whereas 
LSB galaxies show k ~ 0.2 (|de Blok et al. 2001 j) . Other 
analyses including the effects of baryons conclude that n ~ 
0.4 (|Oh et al. 2010bp . This shows that the problem of the 
dark matter central distribution in galaxies is still under 
debate, and there is no clear answer at the moment. Never- 
theless, in this paper we look into an even smaller scale of 
the dark halo, that is, we study the dark matter distribu- 
tion in the vicinity of a supermassive non-rotating black hole 
in order to understand whether or not a particular density 
profile like those mentioned above is favored over the others. 

Supermassive black holes in the center of galaxies con- 
stitute another brick of the problem, because how they were 
formed and how they are fed remains a mayor quest. Assum- 
ing such black holes exist, it is worth asking their influence in 
the distribution of matter around them due to their strong 
gravitational field and determine how they affect in particu- 
lar the dark matter density profile. As described above, the 
distribution of dark matter in galaxies is constructed based 
on simulations obtained from structure formation models, 
and the limitation of such method is the resolution allowed 
by currently used algorithms and implementations, which is 
the reason why the spatial scale of resolution is of the order 
of kpc. We thus explore how a black hole affects the dark 
matter around it and determine whether or not the matter 
distribution is consistent or not with predictions obtained 
from galactic scale analyses and observations. This approach 
conversely, faces the problem of being a rather small scale 
compared to that of dark matter density profiles of galac- 
tic scale, however when studying the problem at black hole 
scale it is possible to infer the behavior of matter in a region 
within lpc, whereas the numerical and observational studies 
above deal with scales of kpc and consider the geometry of 
the center of the galaxy smooth. However we expect that 
our results extend up to the kpc scale. 

An important feature of studying dark matter around a 
black hole is that relativistic fluid dynamics equations on a 
non-flat background space-time are to be solved, unlike the 
most common approach of supermassive black hole growth 
and accretion of dark matter models involving a phase space 
analysis of a Newtonian version of black holes based on the 
seminal accretion models ( |e.g. Li ghtman & Shapir o" 1977| 
IZhao e t al. 2002}. Alternatively, in the relativistic regime 
important bounds on the collisonless dark matter central 
density have been proposed that allow the observed masses 
of central black holes ((Hernandez & Lee 2010]) . In this pa- 
per on the other hand, we model the dark matter as 
a relativistic collisional ideal gas. Even though collisional 
dark matter has been considered to have astrophysical im- 
portance, for instance in non standard models of super- 
massive black hole growth (|Ostriker 20001 IHu et al. 2006]) . 
as a possible solution to the galactic cusp-core problem 
( |Sperg cl & Stci nhardt 2000| and th e corresponding gravi- 
tational collapse (|Moore et al. 2000p . we also focus on the 
limit of a collisionless gas, which is the case consistent with 
the CDM standard approach. 

Two important condition of our model are in turn, 
one is that we do not consider the effects of baryons and 



our model consists on the dark matter component only; 
the second is that we assume the black hole space-time is 
fixed and the dark matter is a test matter field, which is 
justified due to the small accretion rate -even in spheri- 
cal symmetry- happening when the dark matter has pres- 
sure, which is of the order of one solar mass per gigayear. 
( |Guzman & Lor a-Clavij o 2011[ ). We thus solve the fully 
time-dependent relativistic Euler equations for the fluid on 
a Schwarzschild black hole space-time, where the test fluid 
plays the role of dark matter, that is, we implicitly assume 
the dark matter behaves as a smooth matter field, unlike 
the most common approach based on n-body or smoothed 
particle hydrodynamics simulations. We start up modeling 
the dark matter as a relativistic ideal gas with pressure, that 
is, we consider the fluid obeys an ideal gas equation of state, 
an approach that has been used to estimate upper bounds 
on accreted mass of collisional dark matter by supermassive 
black holes ( |Guzman fc LoraTC lavijo 2 011 [ I. We have found 
that there are late-time attractor type of solutions for the 
fluid density profile with constant accretion rate for the cases 
of non-zero pressure. Such late-time scenarios can be found 
starting with rather out of equilibrium initial conditions for 
the fluid in terms of initial inward velocity, initial rest mass 
density and adiabatic index. In this paper we analyze the 
late-time state of the rest mass density profile of such states 
in order to envisage a possible consistent dark matter distri- 
bution at very local scales near a supermassive black hole. 

The paper is organized as follows, in the following sec- 
tion we describe the formulation of Euler's equations we use 
to model the accretion of dark matter, in section [3] we de- 
scribe our numerical methods and in section [4] we present 
our results. Finally in section [5] we discuss our results and 
draw some conclusions. 

2 RELATIVISTIC EULER EQUATIONS 

The system of equations we solve are those of a perfect fluid 
evolving on a fixed curved background space-time, that is 

V„(TH = 0, (1) 
V M (pu") = 0, (2) 

where p is the proper rest mass density, u M is the 4- velocity 
of the fluid and V M is the covariant derivative consistent 
with the 4-metric of the space-time (Misner et al. 1973jl . 
The stress energy tensor of the fluid in these equations is 
given by 

T" v = phu"u + P g' lv . (3) 

where p, h and g^ v are respectively the pressure, the specific 
enthalpy and the 4-metric tensor of the space-time, respec- 
tively. In order to write down explicitly our equations we 
write down the Schwarzschild metric in the 3+1 fashion for 
the coordinate system {t,x % ), that is (jMisner et al. 1973|) : 

da 2 = -(a - p i p i )dt 2 + 2/3 i dx i dt + -y lj dx t dx j , (4) 

where a is the lapse function and /3 8 is the shift vector 
and "fij is the inherited 3-metric on the space-like hyper- 
surfaces. In order to allow the fluid to enter the black 
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Figure 1. (a) We show snapshots of the density from t = to t = 50000M every 250M ; earlier times correspond to the lines in the bottom 
and later times correspond to lines in the upper part approaching the bold black line. The bold black line is actually the superposition 
of various of the late-time snapshots which shows a nearly stationary state. The density profile shows two clear regimes: one near the 
black hole up to r ~ 50M which shows a clear polynomial shape and a second one after r ~ 500M which is nearly constant; these are 
the two domains we use to fit the density profile of the late-time solution, (b) We show a zoom in of a subset of (a) until f = 6000M 
and show that a front shock moving to the right leaves behind the stationary density profile, (c) The velocity of the fluid also reaches 
a steady state with one quarter the speed of light near the horizon, (d) We also show the accretion mass rate measured at the horizon; 
this shows how a steady state is reached at late-times. The initial data correspond to po = 10 — 8 , Vq = —0.1, T = 1.1 and t = 0.8. For 
the accretion mass rate we use the accretion mass rate formula for the spherical accretion case M a cc = —Anr^pW ( v r — — J. 



hole's horizon we use horizon penetrating coordinates with 
Eddington-Finkelstein slices, for which a = , 1 B l = 



\ 2M ( 1 
[ r \l + 2M. 



0,0 and jij 



= diag[l + ^,r 2 , 



- 2 sin 2 6>]. 



For this coordinate choice, equations (II 121) for Q and (j4]) 
written in a balance law type based on the Valencia for- 
mulation of relativistic hydrodynamics ([Ma rti ct al. 19911 
|Banyouls et al. 1997| | Font et al. 2000[> . which for the spher- 
ically symmetric case reads | | Guzman fc Lora-Clavijo 2011[ ) 



du dF r (u) 
dt + dr 



(5) 



where 



D 

Jr 
T 



^/jphW 2 v r 
_ y^fiphW 2 -p- pW) 

§1\ D 



(6) 



In these expressions the 3-velocity measured by an eulerian 
observer v l is denned in terms of the spatial part of the 4- 
velocity of the fluid by v l = |- + ^- , where W is the Lorentz 
factor W — , 1 . . . We close the system of equations 

V 1 1ij v%v3 

with an equation of state p — p(p, e), where e is the specific 
internal energy, defined by h = 1 + e+p/p. For this we choose 
the equation of state of a relativistic ideal gas p = (T— l)pe, 
where F is the adiabatic index. 
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Figure 2. As an illustration of the accuracy of our calculations 
we show the order of self-convergence of p, for the case of initial 
density p = 10~ 8 with T = 1.06 and = -0.1 for e = 0.5. We 
calculate the convergence using the L\ norm of the differences 
between the value of the density for four different resolutions 
Aii = 0.4, Ax 2 = Axi/2, Ax 3 = Ai 2 /2 and Ax 4 = AX3/2, 
corresponding to the numerical calculation of density pi, p%, ps 
and P4 respectively. Then, the self-convergece factor, Q, for these 
three resolutions is calculated from the following expresion as 
2 Qsa = AE1/AE2 where AE± = L 1 {p 1 -p 2 ), AE 2 = L 1 (p 2 -p3) 
and AE3 = Li (p3 — pi). The dominant error is that of the recon- 
struction of variables at the interface cells, which is a piece wise 
constant and therefore first order accurate. We show in this plot 
that Q = 1, that is the relationship AiSj = 2 1 AEi+\ holds. 



2.1 Initial conditions 

In order to solve our system of equations as an initial value 
problem we have to set initial conditions that we choose to 
correspond to out of equilibrium configurations in order to 
let the system of equations to drive the fluid (or not) into 
a relaxed state. Since the nature of dark matter is still un- 
known we have to be as open as possible and choose out of 
equilibrium initial conditions with most free parameters as 
possible. In the first place we choose an initial constant den- 
sity profile with amplitude po = p(t = 0) which we found to 
provide quite out of equilibrium initial conditions. Secondly, 
no initial velocity Vq of the fluid can be assumed, and since 
we start up with out of equilibrium conditions we choose 
this parameter to be rather arbitrary. Also, there is no pre- 
scription -as far as we can tell- about the internal energy of 
a gas model for dark matter, thus we use various values of 
the initial internal energy of the gas eo = t(t = 0). Finally 
the adiabatic index F determines how collisional the gas is, 
and we look for attractor solutions for various values of this 
parameter as well. 



3 NUMERICAL METHODS 

We solve our equations on a spatial domain ranging from 
M < r ^ 10000M, that is, from inside the black hole's 
horizon at r = M where we apply the excision technique 
(Scidel & Suen 1992) up to an external boundary located 
far away from the horizon; the outer boundary in physical 
units for a black mass of M = 1O 9 M0 corresponds to 0.5pc, 
which is a rather smaller scale compared to previous analyses 
focused on dark halo profiles based on structure formation; 



in all our runs we make sure that the external boundary 
is causally disconnected, so that we avoid the potential in- 
fluence of boundary conditions applied in such boundary on 
our results. In order to solve these equations we evolve initial 
data for the three conservative variables using a finite vol- 
ume cell centered discretization of the domain, second order 
accurate High Resolution Shock Capturing method based 
on the HLLE approximate Riemann solver to calculate the 
numerical fluxes at the cell interfaces (|Harten et al. 1983[) 
and a constant piecewise cell reconstruction of variables 
which is first order accurate QLeVeque 1992| ). The evolution 
is worked out using the method of lines using a third or- 
der Runge-Kutta integrator along the time direction. We 
also implemented a density atmosphere in order to avoid 
the divergence of our equations such that the rest mass den- 
sity is always assumed to be p = max(p, p a tm), where we 
choose values between p a tm = 10~ s and p a tm = 10 -12 with 
which we obtain the same physical results and convergence 
of our implementation, and in fact we also have used such 
values as the initial density in some of our numerical ex- 
periments. Furthermore, we have practiced self-convergence 
tests that show our implementation converges to first or- 
der in the primitive and conservative variables as expected 
according to our cell reconstruction algorithm. 



4 RESULTS 

A pressurless gas system corresponds to T = 1, which 
lacks of spherically symmetric late-time stationary solutions 
and is the reason why we do not study such system here 
( |Guzman fc Lora-Clavijo 201 Instead, for T > 1 the fluid 
stabilizes around a situation of constant accretion mass rate 
and in any case the pressureless case can be studied as a 
limit when V — > 1. Therefore, we now focus on this type 
of configurations and analyze their late-time behavior. As 
an example of how a late-time solution looks like we show 
in Fig. [T] snapshots of the rest mass density, velocity of the 
fluid and mass accretion rate for a particular case. In Fig. 
[2] we illustrate the self-convergence of the density for one of 
our runs, which indicates that our implementation is con- 
sistent with the numerical methods and that our results are 
reliable. 

We consider different values of the initial rest mass den- 
sity at initial time po, the adiabatic index V and the internal 
energy at initial time eo. Once the system relaxes, that is, 
once it approaches a stationary type of state with constant 
mass accretion rate and time independent density, pressure 
and velocity, we fit the resulting density profile with a fit- 
ting function of the type A/r K , where k is the parameter 
that determines the slope of the density profile of the late- 
time attractor type solution. Results shown in Fig. [1] suggest 
two different regions for the fitting of the density profile: re- 
gion (I) ranging from the event horizon up to ~ 50Af where 
a power law is evident, and region (II), which ranges from 
approximately 800Af outwards, which within our domain is 
of the order of tenths of parsec. Region I is related to the 
strong field region whereas region II is the one that would 
eventually match the kpc scale and therefore the results to 
be compared with the galactic scale results. These regions 
vary from one simulation to another because the range of 
the parameter space is wide, however we found 50M to be 
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Figure 3. We show the values of re for a given eo = 0.2 and various 
values pa and the adiabatic index V for region I and various very 
different initial conditions. These results show the independence 
of re on po for a given combination of eo and T. We have made sure 
that this happens in all our simulations. We also show the extrap- 
olation corresponding to the dark matter density in geometrical 
units po = 10 -23 which is equivalent to po ~ IOOMq/jjc 3 , for 
a black hole mass M = 10 9 Mq. The transparent points corre- 
spond to the results of our numerical experiments whereas the 
filled symbols correspond to extrapolations of our results to the 
dark matter density. Similar results are obtained for region II. 



an upper bound of region I in all our cases and the range 
800 - 1500M the boundary of region II. 

Our first result is that for given values of V and eo for 
different values of po, the parameter k results to be the same 
as shown in Fig.[3]for region I. That is, the profile exponent 
k does not depend on the initial value of the constant den- 
sity profile po used to start up our simulations. This allows 
one to consider only the parameter space to be the F — eo 
plane and choose a given value of po that provides better 
accuracy and use the resulting slope for other values of such 
initial parameter. In particular we are interested in explor- 
ing the behavior of po that corresponds to a reasonable value 
of the dark matter density. We illustrate this extrapolation 
also in Fig. [3] We have verified that our stationary solu- 
tions reproduce those exact solutions first found by Michel 
(|Michel 1972|) on the accretion of an ideal gas with pressure. 
As far as we can tell, no definitive results on the attractor 
properties of this solution have been discussed before and 
it is interesting that we have found them as the result of 
the evolution of rather general initial conditions. There are 
then two reasons to support our extrapolation: 1) our so- 
lutions show a universal behavior for different values of the 
initial density for various values of V and eo, and 2) Michel's 
solution would be valid for arbitrary values of the density, 
however not verified at such ultralow densities as an attrac- 
tor solution. 

The exploration of the resulting parameter space in- 
dicates that k depends monotonically on F and eo. These 
results are summarized in Fig. [4] for region I and in Fig. [5] 
for region II. For each value of eo we find different values 
of k depending on the value of T. If our fluid represents 
the dark matter, this implies that the preferred slope of the 
late-time solution for a collisional dark matter candidate de- 
pends both, on the initial internal energy and the adiabatic 
index only. Our results indicate that in order to have a less 



Figure 4. We show the values of re in terms of the adiabatic index 
and initial internal energy for region I. The exponent strongly 
depends on the initial conditions. Nevertheless, in the limit of 
r approaching one, that is, the collisionless case, the exponent 
approaches 1.5 independently of eo- 
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Figure 5. The values of re in terms of the adiabatic index and 
initial internal energy are presented. The domain where the den- 
sity is fitted correspond to 800 ^ r/M ^ 1500. As we can see, the 
bigger the values of T the smaller the values of re. 



cuspy profile, that is, smaller k, a higher adiabatic index is 
needed, meaning that the dark matter is more collisional. 

Concerning region I, the extrapolation of our results in 
the limit F — > 1, that is, the collisionless case, prefers an 
approximate exponent k ~ 1.5 independently of the initial 
internal energy as seen in Fig. [4] It is also noticed that for 
bigger values of V smaller values of « are found. 

On the other hand, for region II the results are less 
coincidental. The first finding is that in the limit F — > 1 
the value of k differs for different initial conditions as shown 
in Fig. [5j and no upper bound can be estimated based on 
our results. Nevertheless, we find that pressure suffices to 
provide very small values of k, in fact smaller than 0.1 for 
r ^ 1.1, which is a result consistent with flat density profiles 
related to cored halos. 
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5 CONCLUSIONS AND DISCUSSION 

The results in this paper are well far from the most gen- 
eral case, which would involve the non-spherical accretion of 
matter, on a black hole with non-zero angular momentum 
and possibly a dynamical geometry ruled by the full non- 
linear Einstein's equations. The approximation of the fluid 
to be a test matter field is appropriate in the present case 
due to the small density values we use and to the quick sta- 
bilization of the accretion rate, although the analysis of the 
rate growth of the event horizon in a full non-linear regime 
would certainly provide more precise results and would al- 
low one to track the growth of the black hole event horizon 
mass. 

On the other hand, as a first approach, the spher- 
ical symmetry is justified and is practical in the sense 
that our results work as bounds of more general cases in- 
volving non-zero angular momentum of matter. At this 
respect there are studies considering non spherical ac- 
cretion of matter, for instance {Read & Gilmor c" 20031 
|Holley -Bockcl man fc Sigurdsson 2006[ l, all of them based on 
the phase space analysis focused on the lose-cone refilling 
ruled by Newtonian gravity and therefore an approximation 
that deals well with non-trivial impact parameters of the 
matter but not with the strong gravitational field. Thus a 
study of matter winds with non-zero angular momentum 
around black holes would tell us more about not only den- 
sity profiles of potentially existing attractor in time config- 
urations, but also on the cross section of black holes that 
would better restrict the accretion rates of dark matter. 

Nonetheless, the present first approach results provide 
some important indications. Under the conditions we used, 
the density profile is cuspy only in the region near the black 
hole, whereas the density stabilizes around a nearly con- 
stant density profile in the far region. This result implies 
for instance that for a black hole of mass 1O 9 M0, a den- 
sity profile is nearly flat starting from distances of already 
1000M ~ 0.05 pc on. 

Our main conclusion is that if the profile we found in 
the far region near ~ O.lpc prevails up to scales further 
than those used here, and this is to be compared to galactic 
scale analyses based on structure formation approaches and 
observations, the matter density shows a non-cuspy profile 
for the parameter space we analyzed here. This means that 
a nearly constant dark matter density profile is consistent 
with the presence of a black hole at the price of adding some 
pressure. Moreover, such profile corresponds to an attractor 
in time type of solution, that is, a profile that is obtained 
under a variety of different initial conditions of the matter 
field. 

It is worth to mention here that we are not considering 
the contribution of baryonic matter, however it is impor- 
tant to emphasize that the effects of baryons on the dark 
matter distributions plays an important role for instance in 
the cored shape of dark halos; some of these effects are: 
feedback {Weinberg fc Katz 2002] IValenzuela et al. 2007)) . 
dynamical friction I|E1-Zant et al. 2001 p . triaxial effects 
{Hayashi et al. 2006| |Hayashi et al. 2007| ) among others, 
which should be considered in further relativistic analyses 
like ours. 
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